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Abstract 

In this study, we consider the use of seismic sensors for footstep localization in indoor environ- 
ments. A popular strategy of localization is to use the measured differences in arrival times of 
source signals at multiple pairs of receivers. In the literature, most algorithms that are based 
on time differences of arrival (TDOA) assume that the propagation velocity is a constant as a 
function of the source position, which is valid for air propagation. However, a bounded medium, 
such as a concrete slab, is usually dispersive and damped. This study demonstrates that under 
such conditions, assimilating the concrete slab to a thin plate and considering a Kelvin- Voigt 
damping model, the perceived propagation velocity decreases when the source-sensor distance 
increases. Based on the sign of the measured TDOA (SO-TDOA), a new localization algorithm 
that is adapted to a damped and dispersive medium is proposed here. A simulation and some 
experimental results are included, to define the performance of this SO-TDOA algorithm. 

Keywords: Footstep localization, seismic sensors, time differences of arrival, damping effect, 
dispersive effect, concrete slab, perceived propagation velocity. 



1. Introduction 

For many applications, it is important to obtain location information about a resident in an 
indoor environment. For example, knowing the position of a resident can facilitate the control 
of the heating and air conditioning systems. Existing solutions, however, are intrusive, and they 
do not respect the private life of the resident (e.g., audio or video monitoring [ 1 ]), or they are 
obliging people to keep sensor on their body all the time (e.g., the magneto-inertial navigation 
technique 13). In this study, we propose a new indoor localization algorithm that is not con- 
strained. This new algorithm is based on seismic signal processing. 

The vibration signature of the human footstep on a floor creates an elastic wave that is induced 
by the walking motions. Our goal is to localize footsteps using seismic sensors that are fixed 
on the floor in the indoor environment. Only a few studies have described seismic methods that 
are applicable to footstep localization in an indoor environment. The present techniques can be 
divided into two groups: 

• Techniques based on seismic -wave structures [3 , 4]: with this type of technique, a footstep 
is modeled as a seismic signal composed of P-waves (longitudinal waves) and S-waves 
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(transversal waves) in a three-dimensional environment. Using this assumption, the direc- 
tion of arrival can be determined from the correlation between the signals recorded by a 
three-axis accelerometer. These techniques are not efficient in indoor environments. In- 
deed, the signals recorded indoors by a sensor is a mix of direct and reflected waves (e.g., 
reflections on the edge of the slab, reflections on the furniture and facilities) in an almost 
two-dimensional environment. The time delay between two paths is very short in an indoor 
environment. The distances are only a few meters and the propagation velocity of seismic 
waves is more than 1000m /s in a concrete medium. In addition, elastic waves propagated 
on the floor depend on many factors; e.g., the footwear of the person, the angle of impact 
excitation, the construction of the floor, and the geometrical walking pattern ||5]|6l|7). 

• Techniques based on range delay estimation [8, 9|: these techniques, such as hyperbolic 
localization II 1011 . are based on time differences of arrival (TDOA) and the propagation 
velocity estimation. The propagation velocity is assumed to be constant and independent 
of the source position. In other words, the time of arrival (TOA) depends linearly on the 
source-sensor distance. 

In what follows, we will first discuss the applicability of localization techniques assuming 
a constant propagation velocity for the problem of footstep localization using seismic sensors. 
Indeed, because the various wave components travel at different propagation velocities, footstep 
signals will vary from one receiver location to another. The detected arrival times and the per- 
ceived propagation velocities will closely depend on the attenuation and the dispersion properties 
of the floor. A theoretical study of elastic-wave propagation based on a simplified bending-wave 
equation will be conducted in section|2] This study will show that the perceived propagation ve- 
locity decreases in a floor assimilated to a thin, damped, and dispersive plate if the source-sensor 
distance increases. Analytical and experimental results will also be presented to reinforce this 
conclusion. Therefore localization techniques based on range delay estimation are inadequate 
for our problem. 

A new localization algorithm will be proposed in section [3] This new algorithm takes into ac- 
count the nonconstant propagation velocity and exploits the property that the order of arrival of 
the signals at the sensors is maintained in the dispersive and damped floor being considered. The 
proposed footstep localization algorithm is based on a study of the sign of the time differences of 
the arrival (SO-TDOA). The development of the proposed algorithm will be followed in section 
[4] where we describe simulation results and analyze the performances of the proposed SO-TDOA 
algorithm, as compared with the hyperbolic algorithm that is based on range estimation. Section 
|5]will describe the field tests and provide some experimental results. 



2. Perceived propagation velocity of the seismic signal of a footstep on a floor 

The floor of an indoor environment will be assimilated to a thin damped isotropic plate IfTTl 
[T2l throughout this study. Considering this assumption, the goal of this section is to define 
the influence of the dispersion and the damping effects on the "perceived propagation velocity" 
estimated by a given measuring strategy. 

Consider a plate of thickness h, of infinite extent in the x, y plane. The governing equation for 
the bending motion of a thin undamped plate is |[T3l [Ml : 

d 2 

ph—ru(x,y,t) + DA 2 u(x,y,t) = f (1) 



where u is the transversal displacement, D - E ^fz^ i s tne bending stiffness, E is the Young's 

modulus, <x is the Poisson ratio, p is the mass density, A = J^r + is the Laplacian, and / 
describes the external forces exerted on the plate. Eq. [T] corresponds to the ordinary fiexural 
wave equation. It is satisfied for a thin plate where its thickness h is less than a sixth of the 
wavelength (h < A/ 6). A correction term can also be added in the case of a thick plate, to 
represent the effects of shear stress (although this is not the case in the present study). 
Internal mechanical damping is taken into account by introducing a viscous friction force. This 
friction force is proportional to the time derivative of the strain. Thus Eq. ([TJ for a damped 
medium is given by the Kelvin- Voigt model |fl5l[T6l : 

d 2 I d \ 

ph — u(x,y,t) + D\\ + &-}A 2 u(x,y,t) = f (2) 



dt 2 v ' \ dt. 

Then, the dispersion relation is deduced: 

- a? + a 2 (l -j&aj)k 4 = (3) 

where r\ = -&a> is the dimensionless loss factor that is characteristic of the damping effect, and a 
is a characteristic of the concrete slab, such that a = [m 2 .* -1 ]. This implies that: 

k(a>)= A p(l-j^r», (4) 
V a 

for a low loss factor {§(l> « 1), 

k(to) - + \i §0 ^ = + jki(u)- (5) 

where and k[ are the real and the imaginary parts of the wave number k, respectively, kj is 
known as the attenuation coefficient of the wave in the propagation direction. So the damping 
induces frequency-dependent attenuation (fc/(<y)). The dispersion (k R (a>)) causes a frequency- 
dependent group velocity propagation c g that is given by: 

dco 

c s = tj- - 2 Va yoj. (6) 
ok R 

If the propagation medium is isotropic, the transversal displacement u(x, y, t) depends only on 
the source-sensor distances d. For an initial wave at position (x, y) = (0, 0), the propagative wave 
u(d, t) at another position of distant of d is given by: 

u(d,t) = — f U{Q,u)e miJJ)d - M] du (7) 

= -L f U(0,a))e- kl ^ )d ej [k "^ )d - w,] du (8) 

where U(0, to) is the spectrum of the emitted wave m(0, f). 
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In what follows, we study the variation of the "perceived propagation velocity" with the 
source-sensor distance. Two approaches will be developed. The first calculates the integral in 
Eq. ([8]) using a discrete sum, and then applies a threshold level to detect the time of arrival, and 
thus calculates the "perceived propagation velocity" for a given distance. The second approach 
uses the stationary phase method to evaluate the envelope of the signal in Eq. ([8]), and then to 
determine the relationship between the "perceived propagation velocity" and the source-sensor 
distance. 

We consider a concrete slab of thickness h = 20cm, Young's modulus E = 24A0~ 9 N/m 2 , 
mass density p = 2500kg/m 3 , and Poisson's ratio cr = 0.2, fl6l . Under such conditions, a is 
about 183m 2 /i. This indicated value is useful as an example, because the mechanical properties 
of a material like concrete are known to depend strongly on their composition and how they are 
made. It is also important to note that the expression in Eq. ([8]) is not valid at short source-sensor 
distances (e.g., ~ Im) considering a plate of thickness h ~ 20cm. Indeed, the approximation 
of a thin plate (h < A/ 6) is not valid at these distances because the signal is dominated by high 
frequency components, which is equivalent to a short wavelength (A < Im). 
In the literature, the loss factor of concrete material (rj - &co) can take values from 1CT 3 to 1CT 2 
in the audio frequency range fl6l . Without the lose of generality, we choose § = 10~ 5 s for 
co < \0 A Hz. The approximation {ftco « 1) is then satisfied for the concrete medium. 
For studying the effects of the damped and dispersive medium, we assume that the source signal 
is a Dirac; i.e., U(0, coq) = 1 for all coq. 

2.1. Perceived propagation velocity - integral approximation 

To simulate the received signal u(d, t) at a distance d from the source, an approximation of 
the infinite integral in Eq. <|8j using a discrete finite sum is proposed, with: 

m t) a y , m )* C0S ( kR (^d-^t) (9) 

where co m = 2nf max . f max is fixed at 10kHz and n = 2048 for this simulation. 
It should be noted that the approximation in Eq. ([9]) is not valid for short source-sensor distances 
(d < 3m), because at these distances the signal is dominated by high frequency components that 
are not considered by the finite sum in Eq. (p). 

Figure [T] shows the simulated signals received at d = 5, 10, 15,20m. The amplitude scale is in 
arbitrary units. The TOA can be detected when the signal exceeds a certain threshold A. Figure 
[2] (left) shows the TOA variation as a function of the source-sensor distances for § = 10 -5 and 
§ = 10~ 4 . The threshold value is arbitrarily fixed. Figure [2] (left) also shows the variation of the 
TOA if the perceived propagation velocity is a constant c = lOOOm/s and c = 5000m/ s. 
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Fi gure 1 : Simulated received signal at d = 5, 10, 15, 20/77. 
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Figure 2: Attenuation and dispersion effects on TOA detection (left) and perceived propagation velocity 
(right). 
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Simulation results show that the TOA (t a ) and the source-sensor distance (d) are not linearly 
dependent. The "perceived propagation velocity" c p is defined by: 



c P (d) = — (10) 

ta(d) 

and it is not a constant as a function of the source-sensor distance. However, the order of arrival 
of the signal at the different sensors is maintained (i.e., the TOA increases when the source-sensor 
distance increases). Figure [2] (right side) shows that the perceived propagation velocity appears 
to actually decrease with respect to the propagation distance d. So, if two sensors are placed such 
that sensor 1 is closer to the source k, we have: 

dk\<d k2 => c{d k i) > c(d kl ) (11) 

where du is the distance between the source k and the sensor i, and c{d k i) is the perceived propa- 
gation velocity at the sensors i, and then we obtain: 

tkiidki) = r < t k2 (d k2 ) = d f - (12) 
c(d k i) c(d k2 ) 

where t k i is the TOA detected at sensor i. The TOA detected at the sensor closest to the source is 
the shortest, i.e. : 

If d k \ < d k2 => t kx < t k2 . (13) 

Although the approximation of Eq. (|9]l allows the demonstration of the behavior of the propaga- 
tion wave in a thin plate, the level of approximation, as well as the nature of the approximation, 
barely allows the relationship between d and the TOA to be extracted. To approximate the ex- 
pression of the perceived propagation velocity as a function of the source-sensor distance, we 
use the stationary-phase approximation method, as in the next paragraph. 

2.2. Perceived propagation velocity - stationary-phase approximation 

The stationary-phase method allows the approximation of the evaluation of Eq. in the 
case of a wave packet that propagates in the medium. This leads to the identification of the 
central frequency of the wave packet as a function of d and t lfT31 . This approximation is more 
accurate at around the maximum of the signal. We can write u(d, f) as: 

u{d,f) = — L J U(0, oS)e- k ' (h,)d e^ k ^ d -^duy = j F(oj)e- j ^da) (14) 

where: 

F(aj) = ffi^-fcfcW ( i5) 
V2^ 

O(w) = -[k R (aj)d - tot] (16) 

The stationary phase method consists of expanding €>(<y) in a Taylor series near the point a>o of 
the stationary phase (i.e <S>'{a>o) = 0), keeping only the first two nonzero terms: 

1 „ , 
O(w) =! <D(w ) ( w o)(w - w r (17) 



and approaching F(oj) by F(a>o), the integral in Eq. ( 14 1 can be approached by: 



In 



j|0"(o>o)| 

By inserting the stationary phase condition (<J>'(dL»o) = 0), we get: 

/ d x2 

The envelope of the wave can then be calculated as: 



(18) 
(19) 



(20) 



A(d, f) = 



U(Q,a>o) e _ kl(0)o)d a 2 ff(0,tto)fl*<H4g 



^jd\k' R (oj )\ 



(21) 



t/(o,^) _d_ -^27?- 

f 3/2 

We want to establish the relation between the perceived propagation velocity and the source- 
sensor distance. The proposed approach consists of studying the evolution of the maximum of 
the envelope in time and distance d from the source position. The maximum of the envelope of 
dA 

the signal satisfies — — = 0, then: 
ad 

& 



1-4- 



T d« = 



32a 2 f 3 

and the maximum of the envelope is located at each time at d, which is given by: 

1 



d = 



8a 2 

IT 



(22) 



(23) 



In other terms, the TOA of the maximum of the envelop at a distance d is: 

1 

& 

The "perceived propagation velocity" can then be calculated as: 



(24) 
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(25) 



Eq. ( |25] > shows that the "perceived propagation velocity" is not a constant as a function of 
the source-sensor distance, as it varies like a constant multiplied by (1/d)^ 3 . Moreover, it shows 
that the "perceived propagation velocity" decreases when the source-sensor distance increases. 



These results reinforces those of section 2. 1 which were obtained by a numerical approximation 
of the integral (|8j. 

Figure [3] shows the simulated amplitude of the envelope A(d, f) given by Eq. ( pT| i, as a function 
of time and source-sensor distance {a = 183 and § = 10~ 5 s). It also shows the movement of 
the maximum of the amplitude (thick line) and the movement of a point defined by a constant 
envelope level {I = 200) (thin line) in time and distance. 



A(d.l) 




Figure 3: Envelope of the signal amplitude as function of time and source-sensor distance. Thick line, 
movement of the maximum of the amplitude in time and distance. Thin line, movement of a point defined by 
a constant envelope level (I = 200) in time and distance. 



Experimentally, for large source-sensor distances, high frequencies are severely damped and 
the signal is dominated by the low frequency components. Detection of the maximum suffers 
therefore from high variance. A threshold-based approach can be numerically solved to give the 
shape of the variation of the perceived propagation velocity as a function of the source-sensor 
distances, using the envelope of the signal. The analytical solution of the threshold-based ap- 
proach cannot be easily determined. 

Figure |4] shows the "perceived propagation velocity", as determined by the numerical and 
analytic solutions of the method, solving |j = 0, and by the numerical solution of the method 
solving A(d, t) = 200. This concludes that the "perceived propagation velocity" depends on the 
source-sensor distance. 
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2.3. Experimental results 

To confirm the theoretical relationship between the perceived propagation velocity and source- 
sensor distance, an experimental approach was considered. For the experimental results, we use 
data recorded during indoor tests. The propagation medium considered is a 20 - cm-thick con- 
crete slab covered by linoleum. The sensors (accelerometers) are deployed in a linear array. To 
characterize the propagation for the medium, we used a reproducible source: a ball was dropped 
from a height of 1.50 m several times near a reference sensor go (Figure|5]l. 
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Figure 5: Attenuation and dispersion effects on the perceived propagation velocity: experimental set-up. 
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Figure [6] shows the theoretical "perceived propagation velocity" given by Eq. ( 25 i for & = 
1.5.1(T 5 s and the estimated velocity at each sensor by: 



t a (dd - t a (d ) 

where d, is the distance between sensor i and sensor 0, and t a (di) is the estimated TOA at sensor 
i, as determined by the threshold level on the signal. The choice of the threshold level depends 
on the measured noise level. 



35DD 
3000 
| 2500 



500 






■ 4 




















\ 
























A. ! 


! 










(¥)"'" 






































































* j 






































i 









01 23456769 10 11 12 

Source-sensor distance d [m] 



Figure 6: Attenuation and dispersion effects for the perceived propagation velocity: simulation and exper- 
imental results. (*) estimated propagation velocity for one ball drop; (o) mean of all estimated propagation 
velocities for each distance; (\) error for the estimated propagation velocity related to a 0.1 ms (2 samples) 
error on the measured TOA. 

Figure [6] shows the similarity between the experimental and theoretical variations of the per- 
ceived propagation velocity with distance. It should be noted that the theoretical result corre- 
sponds to the movement of the maximum of the amplitude, and that the experimental result is 
obtained from the movement of the beginning of the signal that exceeds the noise. 
The experimental results confirm again that the "perceived propagation velocity" decreases when 
the source-sensor distance increases in a damped and dispersive thin plate. 

2.4. Conclusion 

To summarize this section, we have shown that the propagation velocity estimation depends 
on the damping and dispersion effects and on the source-sensor distance, and we have shown 
the relation c p (d). Consequently, source localization techniques based on different range esti- 
mations are not applicable. However, we observed that the order of the arrival at the sensors 
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is maintained even in the presence of damping and dispersion. Experimental tests in an indoor 
environment confirmed these results. However, in some cases, and especially when the floor was 
not orthotropic due to the presence of beams in its construction, the order might not be main- 
tained. 

Thus a new algorithm based on the sign of time delay promises good localisation results. In the 
next section we propose a new SO-TDOA algorithm. 



3. Novel SO-TDOA algorithm 

Assuming a damped and dispersive floor, the problem of footstep localization in indoor en- 
vironments cannot be solved using traditional source-localization algorithms based on range es- 
timations, because the perceived propagation velocity depends on the source-sensor distances 
(section [2}. However, the ordering of the arrival time of the signals at different sensors is main- 
tained even in the presence of dissipation and dispersion. In other words, for all of the source 
positions p k and sensor pair (i, j): 



sign (t ki - t k /j = sign 



dki d kj \ , v 
= sign (dki -dkj), 

Cki c kj v 'I 



(27) 



where sign defines the sign operator, and c4,- is the distance between the point pk and the sensor i, 
is the TOA of the signal to the sensor i, and is the perceived propagation velocity at sensor 



i. Eq. ( 27 1 shows that the sign of the time delay is independent of the elastic wave propagation 
velocity in the medium. 

Considering a pair of sensors (/, f) and a point pk, the set S'/ of points that satisfy for all p', € S'/, 



sign [d ki - d k j) = sign (d' u - d' k] ) , 



(28) 



where dki (resp. d',) is the distance between sensor i, and pk (resp. p' k ) is the half space delimited 
by the perpendicular bisectors of the line segment joining the sensors (/, j) and containing pk (see 
Figure |7J. 

Considering now sensors placed in a bounded environment E c ^R 2 . Each sensor is located 
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Fi gure 7; Region separation. 



at a known position gj. The environment is partitioned into Q disjointed regions Rk- Each region 
is limited by the perpendicular bisectors of the line segments joining a pair of sensors. Figure [8] 
illustrates an example of the configuration using 5 sensors in a square room and in a rectangular 
room. 
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From Eq. (27 i and Eq. (28 1, we can deduce the following property. For all points p k and p' k in a 
region and for all pairs of sensors (i, j), we can write: 



VPk, p' k eR k, sign (t ki - f y ) = sign (t' ki - f^) (29) 

The SO-TDOA algorithm consists on region localization. In what follows, we propose to 
characterize each region formed by perpendicular bisector of pairs of sensors. So we will deter- 
mine the number Q of the obtained regions, the coordinates of their centroid point pi, and their 
characteristic vector Zk, as defined below. 

3.1. Region characteristic vector 

Considering all of the sensor pairs (i, j), we can define a characteristic vector z k for each 
region R^ as: 

z k (l) = sign (d u - d kj ) , I = + i ( 30) 

V(i, j) e {(1,2), (1, 3), (2, 3), (1,4), (2,4), ...,(#- UN)}- 

The vector is formed by N(N - l)/2 elements taking values in {+1, 1}. 

Example: Considering the previous example of configuration, the region Ri can be defined by 
the vector z\ of 10 elements, as in Figure[8] 




Figure 8: Example of region separation: 5 sensors in a square room (right) and in a rectangular 
room (left). In the square room, the region R\ is delimited by the perpendicular bisectors of 
the line segments joining pairs (1,5) and (2,4), and the boundary of the environment. z\ is the 
characteristic vector of the region R\ (shaded). 

Remarks: If N sensors are placed in an unbounded plane such that there are no parallel per- 
pendicular bisectors, the number of perpendicular bisectors is equal to N(N - l)/2 (i.e., to the 
number of pairs of sensors). Or considering o nonparallel lines in an unbounded plane, these 
form o(o + l)/2 + 1 regions. The number of regions formed by sensors in an infinite space is 
calculated for o = N(N - l)/2 and is obtained as M R = (N 4 - 2N 3 + 3N 2 - 2A0/8 + 1. M R is 
the upper boundary of the number of regions in a bounded plane. Indeed, the number of regions 
in a bounded plane depends on the number of sensors, their locations, and the room geometry. 
Therefore, there is no simple expression that gives the number of regions formed in a bounded 
plane according to a given sensor configuration. For example, the upper bounds of the number 
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of regions formed with N — 5 sensors is M R = 56. However the number of regions formed in a 
square room is Q = 16, and in a rectangular room, Q = 20 (see Figure [8]). 
The value of ik is in {+1, -\} N ( N - l V 2 if it is considered that the estimated SO-TDOA might be er- 
roneous for some pairs of sensors. As 2 N(N ^^ 2 » Q, only a few values of Zk actually correspond 
to one of the acceptable Q regions. For example, in FigureJH] the sensor pairs (3,4) and (1,2) 
share the same perpendicular bisector, and so the corresponding elements in the characteristic 
vector must have the same value +1 or -1. However, under experimental conditions and with the 
presence of TDOA estimation errors, nonrealistic characteristic vectors can be obtained. Thus, 
using redundancy in the characteristic vector might lead to improved localization performances. 
To estimate the set of regions M r that correspond to a measured characteristic vector z s , we 
choose to minimize the Hamming distance between the measured vector and all of the accept- 
able characteristic vectors Zk- 

M r = arg min d H (z s , z k ), M r c [I..Q] (31) 

ke[l..Q] 

where dn is the Hamming distance measuring the number of components that are different in 
two vectors, 

N(N-V)(2 

d H (z s ,z k )= fe(0©Z*(0), ( 32 ) 

where © is the exclusive or operator (a © b — lifa^Z? else 0, Va,£> € {+1, -1}). The number of 
regions that minimize the Hamming distances to the measured vector can be > 1 in some cases 
where the measured vector z s does not correspond to a realistic region according to the sensor 
configuration considered. This might frequently occur in the presence of TDOA estimation er- 
rors. We denote \M r \ as the cardinal numbers of the set M r . Then to have \M r \ > 1 is possible. 
Below is an example for the region configuration described in Figure [8] 

z 2 = [-1, -1,-1,-1,-1, +1,-1,+1,+1,+1], 

zi = [-1,-1, -1,-1, +1, +1,-1, +1,+1,+1], (33) 

d H (Zz,Zi) =0 + + + 0+1+0 + + + + 0=1. 

Note that two neighboring regions will be "Hamming"-separated by 1 . 



3.2. Region center coordinates 

All points located in the same region Rk are characterized by the same vector z k , as defined 
by Eq. (30i. All of these points will be associated to their centroid pi. Generally, the geometry 
of the sensor location (which can be arbitrary) does not allow a simple analytical calculation 
of the centroid region coordinates to be obtained. We propose to associate each region with its 
centroid, and to develop a simple computer-based approach to determine its coordinates. This 
consists of sampling the space with regular points for location p e ; see Figure [9] For each point 
p e , we compute the characteristic vector z e . Then, all of these points are classified into groups 
by their characterizing vectors. The number of groups obtained is equal to the number Q of 
the total regions formed. The centroid coordinates of one region are obtained by averaging the 
coordinates of all of the points in the same region. This step of the calculation is performed 
only once, when the sensor configuration is fixed. This information is stored and used later to 
determine the source position. 
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Figure 9: Example: region centroid determination. 



3.3. SO-TDOA localization algorithm 

A human footstep generates a seismic signal that is collected at each sensor in the room. To 
localize this footstep, the SO-TDOA algorithm is proposed. It consists of the following steps: 

1. The time of arrival ?„■ of the seismic signal at each sensor i is estimated by a simple thresh- 
old method. This is determined with respect to a common arbitrary time origin IfTTll . 

2. Then the characteristic vector of the source is determined, such that: 

Z s (l) = sign [t si - i sj ) , (34) 

as arranged in Eq. [30] 

3. The set of regions M r c [1..Q] that minimizes the Hamming distance is estimated: 

N(N-l)/2 

M r = arg min V (z s (i) Z k (i)) ■ (35) 

ke[\..M] 

z'=l 

where the cardinal numbers of M r can be more than one (\M,\ > 1). 

4. Finally, the source position localization p s is estimated by: 

where p° r is the centroid of the region r. The source position estimate corresponds to the 
average of the centroids of all of the regions that minimize the Hamming distance to the 
measured vector. This estimator is a heuristic estimator that will be validated in this study 
by simulation results. This point will be investigated in more detail in future studies. 



4. Performances studies 

In this section, we propose to illustrate the robustness of the proposed SO-TDOA algorithm. 
We compare it with the classical hyperbolic localization algorithm. When the perceived propaga- 
tion velocity is assumed to be a constant, the hyperbolic algorithm is one of the best localization 
algorithms. Theoretically, the perceived propagation velocity depends on the source-sensor dis- 
tance in a damped and dispersive medium. We indicated that the order of arrival of the signal 
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is maintained, but we have no access to the value of the propagation velocity in each point of 
a room because it depends on both the attenuation and the dispersion. The values obtained for 
the estimated propagation velocity might be highly variable, especially in the presence of strong 
attenuation. The shape of the variation of the perceived propagation velocity versus distances 



can be as illustrated in Figure 1 1 



In this simulation, we study the performances of the two algorithms when the perceived propaga- 
tion velocity varies, as shown in Figure 11 In concrete, the propagation velocity c can vary from 
some hundred to some thousand meters per second, depending on the mechanical and physical 
properties of the medium fT6l . 

For the presented simulations that are based on the hyperbolic algorithm, c was set in the range of 
500m/ s - 2000m/ s. These values correspond to reasonable experimentally encountered values. 



The simulation steps are given in the next paragraph and summarized in Figure 10 
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Figure 10: Simulation steps. 
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4.1. Simulation steps 



Figure 10 illustrates the different steps of the simulation that was conducted to compare the 
proposed SO-TDOA algorithm with the hyperbolic algorithm lfT0l l9l. N sensors are placed in a 
L x x L y rectangular room, with coordinates gj = [jc,- y,], for all 1 < i < N. 



Inputs 

1. A source position is fixed at p s = [x y], such that < x < L x and < y < L y . 

2. The distances from the source to all of the sensors are calculated, as d S i = \\p s - gi\\z, for 
all 1 < i < N. We note d s , the vector of range differences, such that d s ([) = (d S i - d S jj is 
arranged like in Eq. p4| , 

3. Assuming that the perceived elastic-wave propagation velocity is as given by Figure [TT] 
we calculate the arrival times at the sensors as t S i = jtJ^, for all 1 < i < N. Under the ex- 
perimental conditions, we do not have access to the variation in the perceived propagation 
velocity versus the source-sensor distance. This last closely depends on the properties of 
the propagation medium. 

4. The times of arrival f SI are embedded in an additive zero-mean Gaussian perturbation, 
with variance cr,. We assume that this is white and is independent of the signal or the f J( . 
Experimentally, we have obtained a TOA detection error usually in the range of [0 - l]ms 
(i.e., for 68% in [0 - l]ms this implies cr, - 0.5ms. An error of 1ms in the TOA detection 
implies an error of 0.5m to 4m in the source-sensor distance estimation (for propagation 
velocity c e [500 - 4000]m/s). 

5. All time delays are determined. Let f be the vector of A^(A^ - l)/2 time delays, such that: 

f (I) = t si - f sj (37) 
arranged as in Eq. ( |34) , The vector f is the input of the localization algorithms. 



Algorithms 

Both the hyperbolic and SO-TDOA algorithms take t as their input. The hyperbolic algo- 
rithm requires multiple operations to invert the problem. We generate a grid of points that are 
uniformly distributed in the room, and we search for the point that minimizes the criteria corre- 
sponding to the hyperbolic algorithm and the point that minimizes the criteria of the SO-TDOA 
algorithm: 

6) We generate a regular grid of points p = [x p y p ] uniformly distributed in the room. 

7) For each point of the grid, we calculate the range differences vector d p , such that d p (l) = 



(d P i — dpj^j arranged as in Eq. ( 34 1 and d P i — \\p — g,[|2, for all 1 < i < N. 



8) We estimate the source position by the new algorithm based on sign of time delay estima- 
tion Psistg, and b y the hyperbolic algorithm p SVtyfer . 
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New SO-TDOA algorithm 

8.1. For all of the points p, we calculate the vector z p , such that 

Z P (l) = sign (d p (l)) , \fle[l---N(N- l)/2]. (38) 

8.2. The source position estimated is then given by 

N(N-l)/2 

Aw.= flr £ min 2 ( 2 (0ez p (/)), (39) 
p i=\ 

where £ = signif). 
Hyperbolic algorithm 

8.1. We calculate d = cf, where c is a mean propagation velocity that can be estimated 
beforehand from a known source location and estimated time delay (8) . 

8.2. The source position estimated is then given by 

Pin*,* =argrmn\\d-d p \\ 2 . (40) 




500 - 



1 2 3 4 5 6 7 6 9 10 11 12 13 14 15 
Source-sensor distance d [in] 

Figure 11: Shape of the perceived propagation velocity variation versus distance. 
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4.2. Simulation results 

We consider 9 sensors placed in a room 10m x 10m. A source positions is chosen arbitrarily 
for this study p s = [1 3]m, as in Figure [12] The grid of points needed for the localization 
algorithms is generated using 25 x 25 regular points. The performance index that we use is the 



root mean squared error (RMSE) between the estimated p and the actual position p s , as Eq. (41 




91 " g" 2 "93 

Fi gure 12: Configuration study: Nine sensors positioned in a 10m X 10m room. Source positions (*). 



Figure 13 shows the RMSE of the estimated position as a function of cr, for three different 



shapes of perceived propagation velocity variation (Figure 1 1 1 at the same source position. cr t 



is the standard deviation of the noise simulating TOA estimation errors. Results are obtained by 
averaging over M c = 500 Monte Carlo runs for all of the investigated scenarios, as for Eq. (41 1. 



RMSE(o-,) 



1 



2Jli»»-A(W 



(41) 



Figure 13 shows that the proposed localization algorithm SO-TDOA can achieve good local- 
ization results (RMSE < 1.5m) even at high TOA estimator error (cr, = 1ms) in a room of 
10m x 10m without the need for propagation velocity estimation. 

The hyperbolic algorithm performance depends on the velocity estimation. For different shapes 
of variation of the perceived propagation velocity, the performances of the hyperbolic algorithm 
are changing. We observe that the SO-TDOA algorithm is more robust versus a changing veloc- 
ity . 

Finally, it is important to note that the proposed SO-TDOA algorithm is more rapid and has a 
lower calculation cost compared to the hyperbolic localization algorithm. 
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Figure 13: Performance study at the source position p s 
versus distance, as in Figure 1 1 
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= [1 3]m. Perceived propagation velocity varys 



5. Experiment results 



5.1. Test set-up 

To validate and assess the performances of the newly developed algorithm, we used data 
recorded during a series of indoor tests. The soil is a 20cm concrete slab covered by a tiled 
floor. As shown in Figure 14 nine seismic sensors where placed in a rectangular array on a room 




Figure 14: Experimental environment: nine sensors in a rectangular room. 

3.6m x 5.4m. Two types of sensors (accelerometers) were used: six piezo-electric ceramics fixed 
on the floor, with a weight of 5 kg, and three Colibry SF3000L fixed with double-faced tape 
lfT8l . The seismic data was acquired, digitized, and relayed to a mobile data-recording station 
(YOKOGAWA |[T9l ). The seismic data were sampled at 20kHz. Seven footsteps were monitored 



in the location giving in Figure 14 



5.2. Example of experimental signals 

An example of a footstep signal and its time frequency representation are given in Figure 



15 The signal considered corresponds to footstep PI measured at sensor CI (source-sensor 
distance, 4.15m). The parameters of the short-term Fourier transform are for a Hamming window 
of length 128 (6.4ms), an overlapping segment length of 126, a fast Fourier transform length of 
128, and a sampling frequency of 20kHz. We observe that the time frequency representation 



of the experimental signal in Figure 15 shows similarity to those of the damping and dispersive 



medium response | Appendix A | . This implies that the assumption (for a damping and dispersive 
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Fi gure 15: Signal and short-term Fourier transform of footstep PI, received at sensor CI. 



floor) that is considered in this study conforms to the experimental results. Figure 16 shows an 
example of TOA detection for a seismic footstep signal. 




65 66 67 



70 TOA 71 72 73 74 75 



Figure 16: Example of time-of-arrival detection. A zoom of the signal of footstep PI received at sensor 
CI. 



5.3. Results 

The position estimation errors of the experiment source are given in Table [T] The source 
position estimation error is around some tens of centimeters in a room of 3.6m by 5.4m. 



Footstep 


PI 


P2 


P3 


P4 


P5 


P6 


P7 


Estimation error [m] 


0.68 


0.42 


0.54 


0.33 


0.36 


0.29 


0.57 



Table 1 : Experiment source position estimation error using the SO-TDOA algorithm 
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6. Summary and conclusions 

In this study, we have proposed a new footstep localization algorithm based on the SO- 
TDOA. The SO-TDOA algorithm is easily implemented, and it does not need elastic wave ve- 
locity estimation. Indeed, we first showed that the elastic wave propagation velocity varies im- 
portantly with the source position in an indoor environment, where the floor can be defined as 
a thin damped and dispersive plate. Using techniques based on range estimation, like a hyper- 
bolic algorithm, it is not sufficient to estimate the footstep position using seismic sensors. The 
proposed SO-TDOA algorithm provides good simulated and experimental results (a position es- 
timation error of only some tens of centimeters). In future studies, we will adapt SO-TDOA to 
the dynamic localization of a person in an indoor environment. 



Appendix A. Time-frequency analysis 

In this section, we will simulate the time-frequency response of a signal that is propagated 
in a dissipative and dispersive media. We consider a slab of thickness 0.2m, Young's modulus 
E = 24.1(T 9 A7m 2 , and mass density p = 2500kg /m 3 lfT6l . The slab is rectangular, and of 
width l x and length l y in the direction of x and y, respectively. A sensor is placed at [x r y r ] 
and a source is placed at [x e y e ]m. Assuming that the edges of the slab are sealed. So reflection 
induces a change in direction of the displacement. The received signal is the results of successive 
reflections at the edges. We can represent these reflections from source "images", as for Figure 



A. 17 The x and y positions of the source images are given by the two sets: 




Figure A. 17: Examples of source "images". 



For x: [x e + 2pl x - x e + 2pl x ] {p e [— oo + oo], integer) 
Fory: [y e +2ql y -y e + 2ql y ] {q e [-00 +oo], integer) 

The first series underwent an even number of bounces. The signals from these sources should be 
multiplied by R pq = 1. The second should be multiplied by R pq = — 1. 



22 



We calculate for p < P max and for q < Q max the distance between the source image (p, q) and the 
sensor, denoted by d pq . We deduce the signal received at the sensor using: 



u(t) = —— ^ R pq u(d pq , t) 



(A.l) 



*P9 PA 



where u(d pq , i) is given for d = d pq (according to Eq. d9 



Hdt t) , £ e- k '('^ d cos (k R (^) d - ^t) 



where a> m = 2nf max . f max is fixed as 10kHz and n - 2024 for this simulation. 



(A.2) 



Considering a rectangular room of width l x — 3.6m and length l y - 5 Am, a source position at 
x e = 0.9m, y e = 1.35m (PI) and a sensor position at x r = 0.2m, y r = 5.2m (C7). The parameters 
for the short-term Fourier transform are a Hamming window of length of 128 (6Ams), which 
overlapps a segment length 126, a fast Fourier transform length 128, and frequency sampling 



A.18 



of 20 kHz. The simulated signal and time frequency spectrogram are given in Figures A. 19 and 
for a loss factor § = 10~ 5 . The simul ated ti me-frequency spectrogram for a loss factor of 
10~ 6 s and § = I0~ 4 s are given in Figure A. 20 As compared to the experimental spectrogram 



signal given in Figure [151 the more similar spectrogram signal is given for & = 10" 3 s. 
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